Traffic Dataset
11/18/2020
Input Data
| DateTime | Junction | Vehicles | ID |
|---|---|---|---|
| 2015-11-01 00:00:00 | 1 | 15 | 20151101001 |
| 2015-11-01 01:00:00 | 1 | 13 | 20151101011 |
| 2015-11-01 02:00:00 | 1 | 10 | 20151101021 |
| 2015-11-01 03:00:00 | 1 | 7 | 20151101031 |
| 2015-11-01 04:00:00 | 1 | 9 | 20151101041 |
| 2015-11-01 05:00:00 | 1 | 6 | 20151101051 |
| 2015-11-01 06:00:00 | 1 | 9 | 20151101061 |
| 2015-11-01 07:00:00 | 1 | 8 | 20151101071 |
| 2015-11-01 08:00:00 | 1 | 11 | 20151101081 |
| 2015-11-01 09:00:00 | 1 | 12 | 20151101091 |
| DateTime | Junction | Vehicles | ID | |
|---|---|---|---|---|
| 48111 | 2017-06-30 14:00:00 | 4 | 10 | 20170630144 |
| 48112 | 2017-06-30 15:00:00 | 4 | 14 | 20170630154 |
| 48113 | 2017-06-30 16:00:00 | 4 | 16 | 20170630164 |
| 48114 | 2017-06-30 17:00:00 | 4 | 16 | 20170630174 |
| 48115 | 2017-06-30 18:00:00 | 4 | 17 | 20170630184 |
| 48116 | 2017-06-30 19:00:00 | 4 | 11 | 20170630194 |
| 48117 | 2017-06-30 20:00:00 | 4 | 30 | 20170630204 |
| 48118 | 2017-06-30 21:00:00 | 4 | 16 | 20170630214 |
| 48119 | 2017-06-30 22:00:00 | 4 | 22 | 20170630224 |
| 48120 | 2017-06-30 23:00:00 | 4 | 12 | 20170630234 |
## nextknitr::kable(rbind(head(cars),tail(cars)))
cat("next")| Date | Hour | junction1 | junction2 | junction3 | |
|---|---|---|---|---|---|
| 1 | 2015-11-01 | 0 | 15 | 6 | 9 |
| 2 | 2015-11-01 | 1 | 13 | 6 | 7 |
| 3 | 2015-11-01 | 2 | 10 | 5 | 5 |
| 4 | 2015-11-01 | 3 | 7 | 6 | 1 |
| 5 | 2015-11-01 | 4 | 9 | 7 | 2 |
| 6 | 2015-11-01 | 5 | 6 | 2 | 2 |
| 14587 | 2017-06-30 | 18 | 95 | 34 | 38 |
| 14588 | 2017-06-30 | 19 | 105 | 34 | 33 |
| 14589 | 2017-06-30 | 20 | 96 | 35 | 31 |
| 14590 | 2017-06-30 | 21 | 90 | 31 | 28 |
| 14591 | 2017-06-30 | 22 | 84 | 29 | 26 |
| 14592 | 2017-06-30 | 23 | 78 | 27 | 39 |
## next# knitr::kable(head(traffic.data))
# knitr::kable(head(traffic.data[traffic.data$Junction == 2,]))
# knitr::kable(head(traffic.data[traffic.data$Junction == 3,]))
knitr::kable(head(traffic.data[traffic.data$Junction == 4,1:3]))
knitr::kable(head(cbind(traffic.data[traffic.data$Junction == 1,1:3][1:10,],
traffic.data[traffic.data$Junction == 2,2:3][1:10,],
traffic.data[traffic.data$Junction == 3,2:3][1:10,])))
hold = cbind(traffic.data[traffic.data$Junction == 1,1:3],
traffic.data[traffic.data$Junction == 2,2:3],
traffic.data[traffic.data$Junction == 3,2:3])
cat("next")| DateTime | Junction | Vehicles | |
|---|---|---|---|
| 43777 | 2017-01-01 00:00:00 | 4 | 3 |
| 43778 | 2017-01-01 01:00:00 | 4 | 1 |
| 43779 | 2017-01-01 02:00:00 | 4 | 4 |
| 43780 | 2017-01-01 03:00:00 | 4 | 4 |
| 43781 | 2017-01-01 04:00:00 | 4 | 2 |
| 43782 | 2017-01-01 05:00:00 | 4 | 1 |
| DateTime | Junction | Vehicles | Junction | Vehicles | Junction | Vehicles |
|---|---|---|---|---|---|---|
| 2015-11-01 00:00:00 | 1 | 15 | 2 | 6 | 3 | 9 |
| 2015-11-01 01:00:00 | 1 | 13 | 2 | 6 | 3 | 7 |
| 2015-11-01 02:00:00 | 1 | 10 | 2 | 5 | 3 | 5 |
| 2015-11-01 03:00:00 | 1 | 7 | 2 | 6 | 3 | 1 |
| 2015-11-01 04:00:00 | 1 | 9 | 2 | 7 | 3 | 2 |
| 2015-11-01 05:00:00 | 1 | 6 | 2 | 2 | 3 | 2 |
## nextcreate daily cars
dates = unique(cars$Date)
cars.day = data.frame(Date = dates, Junction1 = rep(NA,length(dates)), Junction2 = rep(NA,length(dates)), Junction3 = rep(NA,length(dates)))
for(i in 1:length(dates)){
cars.day$Junction1[i] = sum(cars$junction1[cars$Date == dates[i]])
cars.day$Junction2[i] = sum(cars$junction2[cars$Date == dates[i]])
cars.day$Junction3[i] = sum(cars$junction3[cars$Date == dates[i]])
}
cars.day$Junction1 = ts(cars.day$Junction1, frequency = 7)
cars.day$Junction2 = ts(cars.day$Junction2, frequency = 7)
cars.day$Junction3 = ts(cars.day$Junction3, frequency = 7)
knitr::kable(rbind(head(cars.day),tail(cars.day)))
plot(cars.day$Junction1, type = "n", ylim = c(0,2250), main = "daily cars at junctions", xlab = "Time (day)", ylab = "cars")
grid()
lines(cars.day$Junction1, col = rainbow(3)[1])
lines(cars.day$Junction2, col = rainbow(3)[2])
lines(cars.day$Junction3, col = rainbow(3)[3])
legend("topleft",legend = c("Junction 1", "Junction 2", "Junction 3"), col = rainbow(3),lty = c(1,1,1))
plot(cars.day$Junction1, type = "n", ylim = c(0,2250), main = "daily cars at junctions 1 & 2", xlab = "Time (day)", ylab = "cars")
grid()
lines(cars.day$Junction1, col = rainbow(3)[1])
lines(cars.day$Junction2, col = rainbow(3)[2])
legend("topleft",legend = c("Junction 1", "Junction 2"), col = rainbow(3)[1:2],lty = c(1,1))| Date | Junction1 | Junction2 | Junction3 | |
|---|---|---|---|---|
| 1 | 2015-11-01 | 327 | 133 | 136 |
| 2 | 2015-11-02 | 546 | 197 | 166 |
| 3 | 2015-11-03 | 544 | 217 | 150 |
| 4 | 2015-11-04 | 498 | 199 | 121 |
| 5 | 2015-11-05 | 464 | 200 | 106 |
| 6 | 2015-11-06 | 446 | 199 | 101 |
| 603 | 2017-06-25 | 1184 | 367 | 335 |
| 604 | 2017-06-26 | 1774 | 587 | 414 |
| 605 | 2017-06-27 | 2187 | 730 | 567 |
| 606 | 2017-06-28 | 2080 | 720 | 503 |
| 607 | 2017-06-29 | 2086 | 696 | 497 |
| 608 | 2017-06-30 | 1883 | 655 | 553 |
acf(cars.day$Junction1, lag.max = 30, main = "Junction 1 ACF")
pacf(cars.day$Junction1, lag.max = 30, main = "Junction 1 PACF")
acf(cars.day$Junction2, lag.max = 30, main = "Junction 2 ACF")
pacf(cars.day$Junction2, lag.max = 30, main = "Junction 2 PACF")
acf(cars.day$Junction3, lag.max = 30, main = "Junction 3 ACF")
pacf(cars.day$Junction3, lag.max = 30, main = "Junction 3 PACF")
create traning data
# remove last month data
nday = nrow(cars.day)
cars.day.train = cars.day[-((nday-29):nday),] # remove last month
cars.day.train$Junction1 = ts(cars.day.train$Junction1,frequency = 7)
cars.day.train$Junction2 = ts(cars.day.train$Junction2,frequency = 7)
cars.day.train$Junction3 = ts(cars.day.train$Junction3,frequency = 7)auto arima
Junction 1
(cars.day.fit1 = auto.arima(cars.day.train$Junction1,
max.p = 5,
max.q = 5,
max.P = 5,
max.Q = 5,
max.order = 10,
max.d = 3,
max.D = 3,
start.p = 1,
start.q = 1,
start.P = 1,
start.Q = 1,
nmodels = 100,
approximation = T))
#quiet(astsa::sarima(cars.day.train$Junction1, p = 1, d = 1, q = 2))## Series: cars.day.train$Junction1
## ARIMA(1,0,1)(0,1,1)[7] with drift
##
## Coefficients:
## ar1 ma1 sma1 drift
## 0.8128 -0.3473 -0.8211 2.1295
## s.e. 0.0386 0.0629 0.0243 0.3114
##
## sigma^2 = 6325: log likelihood = -3310.96
## AIC=6631.91 AICc=6632.02 BIC=6653.65Junction 2
(cars.day.fit2 = auto.arima(cars.day.train$Junction2,
max.p = 3,
max.q = 3,
max.P = 3,
max.Q = 3,
max.order = 5,
max.d = 1,
max.D = 1,
start.p = 1,
start.q = 1,
start.P = 1,
start.Q = 1,
nmodels = 100,
approximation = T))
quiet(astsa::sarima(cars.day.train$Junction2, p = 3, d = 1, q = 3))## Series: cars.day.train$Junction2
## ARIMA(2,0,1)(0,1,1)[7] with drift
##
## Coefficients:
## ar1 ar2 ma1 sma1 drift
## 1.4178 -0.4346 -0.8493 -0.7995 0.7271
## s.e. 0.0652 0.0593 0.0444 0.0256 0.2996
##
## sigma^2 = 778.1: log likelihood = -2711.57
## AIC=5435.15 AICc=5435.3 BIC=5461.23
Junction 3
(cars.day.fit3 = auto.arima(cars.day.train$Junction3,
max.p = 3,
max.q = 3,
max.P = 3,
max.Q = 3,
max.order = 5,
max.d = 1,
max.D = 1,
start.p = 2,
start.q = 2,
start.P = 2,
start.Q = 2,
nmodels = 50,
approximation = T))
quiet(astsa::sarima(cars.day.train$Junction1, p = 1, d = 1, q = 1))## Series: cars.day.train$Junction3
## ARIMA(1,1,1)(0,0,1)[7]
##
## Coefficients:
## ar1 ma1 sma1
## 0.6969 -0.9755 0.0675
## s.e. 0.0354 0.0125 0.0420
##
## sigma^2 = 8704: log likelihood = -3435.03
## AIC=6878.07 AICc=6878.14 BIC=6895.5
find seasonal model
Junction 1
cat("model 1")
(model1.1 = quiet(astsa::sarima(cars.day.train$Junction1, p = 1, d = 1, q = 2, P = 1, D = 0, Q = 0, S = 7)))
#AIC 6953.23
# model1.forecast2 = forecast(model1.1, h = 30)
# rmse1.1 = rmse(model1.forecast1$mean, data.1.full$J1[577:606])
# 114.3474
cat("model 2")
(model1.2 = quiet(astsa::sarima(cars.day.train$Junction1, p = 1, d = 1, q = 2, P = 0, D = 0, Q = 1, S = 7)))
#AIC 7297.4
cat("model 3")
(model1.3 = quiet(astsa::sarima(cars.day.train$Junction1, p = 1, d= 0, q = 1, P = 1, D = 1, Q = 0, S = 7)))
#AIC 6744.5
cat("model 4")
(model1.4 = quiet(astsa::sarima(cars.day.train$Junction1, p = 1, d= 1, q = 0, P = 2, D = 1, Q = 0, S = 7)))
#AIC 6749.55
cat("model 5")
(model1.5 = quiet(astsa::sarima(cars.day.train$Junction1, p = 1, d= 1, q = 0, P = 2, D = 0, Q = 0, S = 7)))
#AIC 6897.51
data.1 = data.frame(J1 = cars.day.train$Junction1[-c(1,2)],
J2L1 = cars.day.train$Junction2[-c(1,578)],
J2L2 = cars.day.train$Junction2[-c(577,578)],
J3L1 = cars.day.train$Junction3[-c(1,578)],
J3L2 = cars.day.train$Junction3[-c(577,578)])
data.1.full = data.frame(J1 = cars.day$Junction1[-c(1,2)],
J2L1 = cars.day$Junction2[-c(1,608)],
J2L2 = cars.day$Junction2[-c(607,608)],
J3L1 = cars.day$Junction3[-c(1,608)],
J3L2 = cars.day$Junction3[-c(607,608)])
cat("model 6")
# sarima(data.1$J1, p = 1, d= 1, q = 0, P = 2, D = 0, Q = 0, S = 7,
# xreg = cbind(data.1$J2L1,data.1$J3L1))
(model1.6 = quiet(astsa::sarima(data.1$J1,
p = 1,
d= 1,
q = 0,
P = 2,
D = 0,
Q = 0,
S = 7,
xreg = data.1[,-1])))
#AIC 6871.87
# model1.forecast1 = forecast(model1.6, h = 30)
# rmse1.1 = rmse(model1.forecast1$mean, data.1.full$J1[577:606])
# 114.3474
cat("model 7")
(model1.7 = quiet(astsa::sarima(data.1$J1,
p = 1,
d= 1,
q = 1,
P = 2,
D = 0,
Q = 0,
S = 7,
xreg = data.1[,-1])))
#AIC 6821.64
cat("model 8")
(model1.8 = quiet(astsa::sarima(data.1$J1,
p = 2,
d= 1,
q = 1,
P = 2,
D = 0,
Q = 0,
S = 7,
xreg = data.1[,-1])))
#AIC 6802.72
cat("model 9")
(model1.9 = quiet(astsa::sarima(data.1$J1,
p = 2,
d= 0,
q = 1,
P = 2,
D = 1,
Q = 0,
S = 7,
xreg = data.1[,-1])))
#AIC 6675.71
cat("model 10")
(model1.10 = quiet(astsa::sarima(data.1$J1,
p = 2,
d= 1,
q = 1,
P = 2,
D = 1,
Q = 0,
S = 7,
xreg = data.1[,-1])))
#AIC 6672.2
cat("model 11")
(model1.11 = quiet(astsa::sarima(data.1$J1,
p = 3,
d= 1,
q = 1,
P = 2,
D = 1,
Q = 0,
S = 7,
xreg = data.1[,-1])))
#AIC 6666.58
cat("model 12")
(model1.12 = quiet(astsa::sarima(data.1$J1,
p = 2,
d= 1,
q = 3,
P = 2,
D = 1,
Q = 0,
S = 7,
xreg = data.1[,-1])))
#AIC 6666.59
cat("model 13\n\n")
(model1.13<-Arima(data.1$J1,
order=c(1, 0, 1),
seasonal=list(order=c(0, 1, 1), period=7),
include.drift=T,
xreg = as.matrix(data.1[,-1])))
model1.forecast1 = forecast(model1.13, h = 30, xreg = as.matrix(data.1.full[577:606,-1]))
rmse1.1 = rmse(model1.forecast1$mean, data.1.full$J1[577:606])
# 109.0436
#AIC 6610.41
cat("model 14n\n")
(model1.14<-Arima(data.1$J1,
order=c(1, 0, 1),
seasonal=list(order=c(0, 1, 1), period=7),
include.drift=T,
xreg = as.matrix(data.1[,-c(1,3,4,5)])))
#AIC 6604.65
#plot(data.1$J2L1, data.1$J1)## model 1$fit
##
## Call:
## arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, Q), period = S),
## xreg = constant, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc,
## REPORT = 1, reltol = tol))
##
## Coefficients:
## ar1 ma1 ma2 sar1 constant
## 0.6072 -1.1318 0.1318 0.8665 2.1810
## s.e. 0.0632 0.0774 0.0773 0.0219 0.3239
##
## sigma^2 estimated as 9633: log likelihood = -3470.62, aic = 6953.23
##
## $degrees_of_freedom
## [1] 572
##
## $ttable
## Estimate SE t.value p.value
## ar1 0.6072 0.0632 9.6104 0.0000
## ma1 -1.1318 0.0774 -14.6256 0.0000
## ma2 0.1318 0.0773 1.7054 0.0887
## sar1 0.8665 0.0219 39.6120 0.0000
## constant 2.1810 0.3239 6.7330 0.0000
##
## $AIC
## [1] 12.05066
##
## $AICc
## [1] 12.05084
##
## $BIC
## [1] 12.09597
##
## model 2$fit
##
## Call:
## arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, Q), period = S),
## xreg = constant, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc,
## REPORT = 1, reltol = tol))
##
## Coefficients:
## ar1 ma1 ma2 sma1 constant
## 0.1312 -0.5222 -0.4778 0.5960 2.1702
## s.e. 0.0651 0.0551 0.0549 0.0294 0.0884
##
## sigma^2 estimated as 17586: log likelihood = -3642.7, aic = 7297.4
##
## $degrees_of_freedom
## [1] 572
##
## $ttable
## Estimate SE t.value p.value
## ar1 0.1312 0.0651 2.0165 0.0442
## ma1 -0.5222 0.0551 -9.4789 0.0000
## ma2 -0.4778 0.0549 -8.6971 0.0000
## sma1 0.5960 0.0294 20.2882 0.0000
## constant 2.1702 0.0884 24.5467 0.0000
##
## $AIC
## [1] 12.64714
##
## $AICc
## [1] 12.64732
##
## $BIC
## [1] 12.69245
##
## model 3$fit
##
## Call:
## arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, Q), period = S),
## xreg = constant, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc,
## REPORT = 1, reltol = tol))
##
## Coefficients:
## ar1 ma1 sar1 constant
## 0.7900 -0.3288 -0.5089 2.2030
## s.e. 0.0393 0.0579 0.0367 1.1115
##
## sigma^2 estimated as 7724: log likelihood = -3367.25, aic = 6744.5
##
## $degrees_of_freedom
## [1] 567
##
## $ttable
## Estimate SE t.value p.value
## ar1 0.7900 0.0393 20.0824 0.000
## ma1 -0.3288 0.0579 -5.6820 0.000
## sar1 -0.5089 0.0367 -13.8661 0.000
## constant 2.2030 1.1115 1.9820 0.048
##
## $AIC
## [1] 11.81174
##
## $AICc
## [1] 11.81186
##
## $BIC
## [1] 11.8498
##
## model 4$fit
##
## Call:
## arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, Q), period = S),
## include.mean = !no.constant, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc,
## REPORT = 1, reltol = tol))
##
## Coefficients:
## ar1 sar1 sar2
## -0.3676 -0.6937 -0.3084
## s.e. 0.0391 0.0400 0.0401
##
## sigma^2 estimated as 7964: log likelihood = -3370.77, aic = 6749.55
##
## $degrees_of_freedom
## [1] 567
##
## $ttable
## Estimate SE t.value p.value
## ar1 -0.3676 0.0391 -9.3900 0
## sar1 -0.6937 0.0400 -17.3483 0
## sar2 -0.3084 0.0401 -7.6968 0
##
## $AIC
## [1] 11.84131
##
## $AICc
## [1] 11.84139
##
## $BIC
## [1] 11.87181
##
## model 5$fit
##
## Call:
## arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D, Q), period = S),
## xreg = constant, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc,
## REPORT = 1, reltol = tol))
##
## Coefficients:
## ar1 sar1 sar2 constant
## -0.3652 0.4476 0.5056 2.4760
## s.e. 0.0391 0.0357 0.0365 45.8665
##
## sigma^2 estimated as 8697: log likelihood = -3443.76, aic = 6897.51
##
## $degrees_of_freedom
## [1] 573
##
## $ttable
## Estimate SE t.value p.value
## ar1 -0.3652 0.0391 -9.3386 0.000
## sar1 0.4476 0.0357 12.5358 0.000
## sar2 0.5056 0.0365 13.8638 0.000
## constant 2.4760 45.8665 0.0540 0.957
##
## $AIC
## [1] 11.95409
##
## $AICc
## [1] 11.95421
##
## $BIC
## [1] 11.99186
##
## model 6$fit
##
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
## Q), period = S), xreg = xreg, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc,
## REPORT = 1, reltol = tol))
##
## Coefficients:
## ar1 sar1 sar2 J2L1 J2L2 J3L1 J3L2
## -0.3823 0.4515 0.4979 0.2207 -0.2148 -0.0083 -0.0106
## s.e. 0.0404 0.0359 0.0368 0.1195 0.1171 0.0331 0.0330
##
## sigma^2 estimated as 8587: log likelihood = -3427.93, aic = 6871.87
##
## $degrees_of_freedom
## [1] 568
##
## $ttable
## Estimate SE t.value p.value
## ar1 -0.3823 0.0404 -9.4714 0.0000
## sar1 0.4515 0.0359 12.5916 0.0000
## sar2 0.4979 0.0368 13.5462 0.0000
## J2L1 0.2207 0.1195 1.8458 0.0654
## J2L2 -0.2148 0.1171 -1.8336 0.0672
## J3L1 -0.0083 0.0331 -0.2500 0.8026
## J3L2 -0.0106 0.0330 -0.3221 0.7475
##
## $AIC
## [1] 11.95108
##
## $AICc
## [1] 11.95142
##
## $BIC
## [1] 12.01166
##
## model 7$fit
##
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
## Q), period = S), xreg = xreg, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc,
## REPORT = 1, reltol = tol))
##
## Coefficients:
## ar1 ma1 sar1 sar2 J2L1 J2L2 J3L1 J3L2
## 0.5188 -0.9898 0.4956 0.4412 0.3452 0.1263 -0.0006 0.027
## s.e. 0.0429 0.0048 0.0375 0.0385 0.1358 0.1304 0.0341 0.034
##
## sigma^2 estimated as 7880: log likelihood = -3401.82, aic = 6821.64
##
## $degrees_of_freedom
## [1] 567
##
## $ttable
## Estimate SE t.value p.value
## ar1 0.5188 0.0429 12.0839 0.0000
## ma1 -0.9898 0.0048 -205.2275 0.0000
## sar1 0.4956 0.0375 13.2063 0.0000
## sar2 0.4412 0.0385 11.4501 0.0000
## J2L1 0.3452 0.1358 2.5421 0.0113
## J2L2 0.1263 0.1304 0.9688 0.3331
## J3L1 -0.0006 0.0341 -0.0168 0.9866
## J3L2 0.0270 0.0340 0.7948 0.4270
##
## $AIC
## [1] 11.86373
##
## $AICc
## [1] 11.86417
##
## $BIC
## [1] 11.93188
##
## model 8$fit
##
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
## Q), period = S), xreg = xreg, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc,
## REPORT = 1, reltol = tol))
##
## Coefficients:
## ar1 ar2 ma1 sar1 sar2 J2L1 J2L2 J3L1 J3L2
## 0.433 0.2043 -0.9928 0.4774 0.4735 0.3400 -0.0389 0.0046 0.0085
## s.e. 0.044 0.0438 0.0038 0.0369 0.0380 0.1271 0.1269 0.0336 0.0336
##
## sigma^2 estimated as 7582: log likelihood = -3391.36, aic = 6802.72
##
## $degrees_of_freedom
## [1] 566
##
## $ttable
## Estimate SE t.value p.value
## ar1 0.4330 0.0440 9.8509 0.0000
## ar2 0.2043 0.0438 4.6667 0.0000
## ma1 -0.9928 0.0038 -259.0924 0.0000
## sar1 0.4774 0.0369 12.9338 0.0000
## sar2 0.4735 0.0380 12.4478 0.0000
## J2L1 0.3400 0.1271 2.6753 0.0077
## J2L2 -0.0389 0.1269 -0.3066 0.7593
## J3L1 0.0046 0.0336 0.1376 0.8906
## J3L2 0.0085 0.0336 0.2519 0.8012
##
## $AIC
## [1] 11.83083
##
## $AICc
## [1] 11.83138
##
## $BIC
## [1] 11.90655
##
## model 9$fit
##
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
## Q), period = S), xreg = xreg, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc,
## REPORT = 1, reltol = tol))
##
## Coefficients:
## ar1 ar2 ma1 sar1 sar2 J2L1 J2L2 J3L1
## 0.8444 0.0069 -0.4156 -0.6663 -0.2948 0.3507 0.0019 -0.0013
## s.e. 0.0947 0.0729 0.0947 0.0416 0.0407 0.1275 0.1260 0.0336
## J3L2
## 0.0004
## s.e. 0.0335
##
## sigma^2 estimated as 6996: log likelihood = -3327.86, aic = 6675.71
##
## $degrees_of_freedom
## [1] 560
##
## $ttable
## Estimate SE t.value p.value
## ar1 0.8444 0.0947 8.9185 0.0000
## ar2 0.0069 0.0729 0.0941 0.9251
## ma1 -0.4156 0.0947 -4.3886 0.0000
## sar1 -0.6663 0.0416 -16.0083 0.0000
## sar2 -0.2948 0.0407 -7.2411 0.0000
## J2L1 0.3507 0.1275 2.7494 0.0062
## J2L2 0.0019 0.1260 0.0153 0.9878
## J3L1 -0.0013 0.0336 -0.0394 0.9686
## J3L2 0.0004 0.0335 0.0110 0.9912
##
## $AIC
## [1] 11.73236
##
## $AICc
## [1] 11.73293
##
## $BIC
## [1] 11.8087
##
## model 10$fit
##
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
## Q), period = S), xreg = xreg, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc,
## REPORT = 1, reltol = tol))
##
## Coefficients:
## ar1 ar2 ma1 sar1 sar2 J2L1 J2L2 J3L1
## 0.4589 0.2249 -1.0000 -0.6481 -0.2838 0.2827 -0.0611 -0.0059
## s.e. 0.0437 0.0431 0.0047 0.0414 0.0409 0.1256 0.1253 0.0337
## J3L2
## -0.0024
## s.e. 0.0337
##
## sigma^2 estimated as 7030: log likelihood = -3326.1, aic = 6672.2
##
## $degrees_of_freedom
## [1] 559
##
## $ttable
## Estimate SE t.value p.value
## ar1 0.4589 0.0437 10.5116 0.0000
## ar2 0.2249 0.0431 5.2188 0.0000
## ma1 -1.0000 0.0047 -211.4468 0.0000
## sar1 -0.6481 0.0414 -15.6632 0.0000
## sar2 -0.2838 0.0409 -6.9340 0.0000
## J2L1 0.2827 0.1256 2.2502 0.0248
## J2L2 -0.0611 0.1253 -0.4881 0.6256
## J3L1 -0.0059 0.0337 -0.1757 0.8606
## J3L2 -0.0024 0.0337 -0.0718 0.9428
##
## $AIC
## [1] 11.74683
##
## $AICc
## [1] 11.74739
##
## $BIC
## [1] 11.82327
##
## model 11$fit
##
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
## Q), period = S), xreg = xreg, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc,
## REPORT = 1, reltol = tol))
##
## Coefficients:
## ar1 ar2 ar3 ma1 sar1 sar2 J2L1 J2L2
## 0.4262 0.1714 0.1190 -1.0000 -0.6629 -0.2953 0.3334 -0.0248
## s.e. 0.0446 0.0471 0.0429 0.0047 0.0414 0.0409 0.1255 0.1246
## J3L1 J3L2
## 0.0007 -0.0023
## s.e. 0.0336 0.0334
##
## sigma^2 estimated as 6938: log likelihood = -3322.29, aic = 6666.58
##
## $degrees_of_freedom
## [1] 558
##
## $ttable
## Estimate SE t.value p.value
## ar1 0.4262 0.0446 9.5476 0.0000
## ar2 0.1714 0.0471 3.6410 0.0003
## ar3 0.1190 0.0429 2.7756 0.0057
## ma1 -1.0000 0.0047 -214.3363 0.0000
## sar1 -0.6629 0.0414 -16.0078 0.0000
## sar2 -0.2953 0.0409 -7.2106 0.0000
## J2L1 0.3334 0.1255 2.6576 0.0081
## J2L2 -0.0248 0.1246 -0.1994 0.8420
## J3L1 0.0007 0.0336 0.0218 0.9826
## J3L2 -0.0023 0.0334 -0.0686 0.9454
##
## $AIC
## [1] 11.73693
##
## $AICc
## [1] 11.73762
##
## $BIC
## [1] 11.82102
##
## model 12$fit
##
## Call:
## stats::arima(x = xdata, order = c(p, d, q), seasonal = list(order = c(P, D,
## Q), period = S), xreg = xreg, transform.pars = trans, fixed = fixed, optim.control = list(trace = trc,
## REPORT = 1, reltol = tol))
##
## Coefficients:
## ar1 ar2 ma1 ma2 ma3 sar1 sar2 J2L1
## -0.1167 0.8114 -0.4463 -0.9757 0.4221 -0.6771 -0.2916 0.3528
## s.e. 0.0621 0.0488 0.0741 0.0457 0.0590 0.0433 0.0412 0.1274
## J2L2 J3L1 J3L2
## -0.0029 -0.0028 0.0047
## s.e. 0.1243 0.0337 0.0335
##
## sigma^2 estimated as 6914: log likelihood = -3321.29, aic = 6666.59
##
## $degrees_of_freedom
## [1] 557
##
## $ttable
## Estimate SE t.value p.value
## ar1 -0.1167 0.0621 -1.8780 0.0609
## ar2 0.8114 0.0488 16.6144 0.0000
## ma1 -0.4463 0.0741 -6.0209 0.0000
## ma2 -0.9757 0.0457 -21.3353 0.0000
## ma3 0.4221 0.0590 7.1540 0.0000
## sar1 -0.6771 0.0433 -15.6260 0.0000
## sar2 -0.2916 0.0412 -7.0873 0.0000
## J2L1 0.3528 0.1274 2.7701 0.0058
## J2L2 -0.0029 0.1243 -0.0233 0.9814
## J3L1 -0.0028 0.0337 -0.0840 0.9331
## J3L2 0.0047 0.0335 0.1391 0.8894
##
## $AIC
## [1] 11.73695
##
## $AICc
## [1] 11.73779
##
## $BIC
## [1] 11.82869
##
## model 13
##
## Series: data.1$J1
## Regression with ARIMA(1,0,1)(0,1,1)[7] errors
##
## Coefficients:
## ar1 ma1 sma1 drift J2L1 J2L2 J3L1 J3L2
## 0.8234 -0.4002 -0.8220 1.9402 0.3404 -0.0569 -0.0028 0.003
## s.e. 0.0390 0.0670 0.0249 0.3212 0.1251 0.1191 0.0354 0.035
##
## sigma^2 = 6300: log likelihood = -3296.21
## AIC=6610.41 AICc=6610.74 BIC=6649.51
## model 14n
## Series: data.1$J1
## Regression with ARIMA(1,0,1)(0,1,1)[7] errors
##
## Coefficients:
## ar1 ma1 sma1 drift xreg
## 0.8220 -0.4012 -0.8206 1.9064 0.3315
## s.e. 0.0394 0.0669 0.0248 0.3123 0.1228
##
## sigma^2 = 6270: log likelihood = -3296.32
## AIC=6604.65 AICc=6604.8 BIC=6630.71
\[ (1-0.822B)(1+0.82B^7)\nabla^{7}J_{1,t} = 1.906 +(1-0.401B)w_{1,t} + J_{2,t-1} \]
J1 model evaluation
hist(model1.14$residuals, main = "histogram of residuals", freq = F, breaks = 35)
lines(-400:400,dnorm(-400:400, mean = mean(model1.14$residuals), sd = sd(model1.14$residuals)), col = "red")
plot(model1.14$residuals, main = "residual plot", type = "l")
adf.test(model1.14$residuals, alternative = "stationary")## Warning in adf.test(model1.14$residuals, alternative = "stationary"): p-value
## smaller than printed p-valueacf(model1.14$residuals, main = "model 14 residual ACF")
pacf(model1.14$residuals, main = "model 14 residual PACF")
model1.forecast1 = forecast(cars.day.fit1, h = 30)
plot(model1.forecast1, xlim = c(75,90))
lines(cars.day$Junction1)
model1.forecast2 = forecast(model1.14, h = 30, xreg = data.1.full$J2L1[577:606])
plot(model1.forecast2, xlim = c(560,610))
lines(data.1.full$J1)
rmse1.1 = rmse(model1.forecast1$mean, data.1.full$J1[577:606])
# 114.3474
rmse1.2 = rmse(model1.forecast2$mean, data.1.full$J1[577:606])
# 109.6458
# mean(data.1.full$J1[577:606])
# 1761.767##
## Augmented Dickey-Fuller Test
##
## data: model1.14$residuals
## Dickey-Fuller = -8.2807, Lag order = 8, p-value = 0.01
## alternative hypothesis: stationary
Junction 2
#quiet(astsa::sarima(cars.day.train$Junction2, p = 3, d = 1, q = 3))
# ARIMA(2,0,1)(0,1,1)[7] with drift
# AIC 5435.15
data.2 = data.frame(J2 = cars.day.train$Junction2[-c(1,2)],
J1L1 = cars.day.train$Junction1[-c(1,578)],
J1L2 = cars.day.train$Junction1[-c(577,578)],
J3L1 = cars.day.train$Junction3[-c(1,578)],
J3L2 = cars.day.train$Junction3[-c(577,578)])
data.2.full = data.frame(J2 = cars.day$Junction2[-c(1,2)],
J1L1 = cars.day$Junction1[-c(1,608)],
J1L2 = cars.day$Junction1[-c(607,608)],
J3L1 = cars.day$Junction3[-c(1,608)],
J3L2 = cars.day$Junction3[-c(607,608)])
# (cars.day.fit2 = auto.arima(data.2$J2,
# max.p = 3,
# max.q = 3,
# max.P = 3,
# max.Q = 3,
# max.order = 5,
# max.d = 1,
# max.D = 1,
# start.p = 2,
# start.q = 1,
# start.P = 1,
# start.Q = 1,
# nmodels = 100,
# approximation = T,
# seasonal = T))
(model2.1 = Arima(data.2$J2,
order=c(2, 0, 1),
seasonal=list(order=c(0, 1, 1), period=7),
include.drift=T,
xreg = as.matrix(data.2[,-c(1)])))
#AIC 5415.31
(model2.2 = Arima(data.2$J2,
order=c(2, 0, 1),
seasonal=list(order=c(0, 1, 1), period=7),
include.drift=T,
xreg = as.matrix(data.2[,-c(1,3)])))
#AIC 5413.33
(model2.3 = Arima(data.2$J2,
order=c(2, 0, 1),
seasonal=list(order=c(0, 1, 1), period=7),
include.drift=T,
xreg = as.matrix(data.2[,-c(1,3,5)])))
#AIC 5411.81## Series: data.2$J2
## Regression with ARIMA(2,0,1)(0,1,1)[7] errors
##
## Coefficients:
## ar1 ar2 ma1 sma1 drift J1L1 J1L2 J3L1 J3L2
## 1.3769 -0.3915 -0.8586 -0.7987 0.6454 0.0301 0.0019 0.0260 0.0080
## s.e. 0.0598 0.0558 0.0365 0.0264 0.3212 0.0157 0.0150 0.0122 0.0122
##
## sigma^2 = 771.7: log likelihood = -2697.66
## AIC=5415.31 AICc=5415.71 BIC=5458.75
## Series: data.2$J2
## Regression with ARIMA(2,0,1)(0,1,1)[7] errors
##
## Coefficients:
## ar1 ar2 ma1 sma1 drift J1L1 J3L1 J3L2
## 1.3771 -0.3917 -0.8580 -0.7990 0.6500 0.0299 0.0260 0.0083
## s.e. 0.0600 0.0560 0.0363 0.0263 0.3193 0.0157 0.0122 0.0119
##
## sigma^2 = 770.3: log likelihood = -2697.66
## AIC=5413.33 AICc=5413.65 BIC=5452.42
## Series: data.2$J2
## Regression with ARIMA(2,0,1)(0,1,1)[7] errors
##
## Coefficients:
## ar1 ar2 ma1 sma1 drift J1L1 J3L1
## 1.3760 -0.3912 -0.8552 -0.7990 0.6554 0.0292 0.0278
## s.e. 0.0611 0.0569 0.0373 0.0263 0.3151 0.0157 0.0119
##
## sigma^2 = 769.6: log likelihood = -2697.9
## AIC=5411.81 AICc=5412.06 BIC=5446.56Junction 2 forecasts
# ARIMA(2,0,1)(0,1,1)[7] with drift
model2.forecast1 = forecast(cars.day.fit2, h = 30)
plot(model2.forecast1, xlim = c(75,90))
lines(cars.day$Junction2)
(rmse2.1 = rmse(model2.forecast1$mean, data.2.full[577:606,1]))
# 43.50
# AIC 5435.15
model2.forecast2 = forecast(model2.1, h = 30,
xreg = as.matrix(data.2.full[577:606,-c(1)]))
plot(model2.forecast2, xlim = c(560,610))
lines(data.2.full$J2)
(rmse2.2 = rmse(model2.forecast2$mean, data.2.full[577:606,1]))
# 43.36
model2.forecast3 = forecast(model2.3, h = 30,
xreg = as.matrix(data.2.full[577:606,-c(1,3,5)]))
plot(model2.forecast3, xlim = c(560,610))
lines(data.2.full$J2)
(rmse2.2 = rmse(model2.forecast3$mean, data.2.full[577:606,1]))
#42.88## [1] 43.50186
## [1] 43.36398
## [1] 42.88388
Junction 3
#quiet(astsa::sarima(cars.day.train$Junction1, p = 1, d = 1, q = 1))
#ARIMA(1,1,1)(0,0,1)[7]
# AIC 6878.07
data.3 = data.frame(J3 = cars.day.train$Junction3[-c(1,2)],
J1L1 = cars.day.train$Junction1[-c(1,578)],
J1L2 = cars.day.train$Junction1[-c(577,578)],
J2L1 = cars.day.train$Junction2[-c(1,578)],
J2L2 = cars.day.train$Junction2[-c(577,578)])
data.3.full = data.frame(J3 = cars.day$Junction3[-c(1,2)],
J1L1 = cars.day$Junction1[-c(1,608)],
J1L2 = cars.day$Junction1[-c(607,608)],
J2L1 = cars.day$Junction2[-c(1,608)],
J2L2 = cars.day$Junction2[-c(607,608)])
# cars.day.fit3
# ARIMA(1,1,1)(0,0,1)[7]
# AIC=6878.07
(model3.1 = Arima(data.3$J3,
order=c(1, 1, 1),
seasonal=list(order=c(0, 0, 1), period=7),
include.drift=T,
xreg = as.matrix(data.3[,-c(1)]),
method = "ML"))
#AIC 6857.59
(model3.2 = Arima(data.3$J3,
order=c(1, 1, 1),
seasonal=list(order=c(0, 0, 1), period=7),
include.drift=T,
xreg = as.matrix(data.3[,-c(1,3,5)]),
method = "ML"))
#AIC 6854.43
(model3.3 = Arima(data.3$J3,
order=c(1, 1, 1),
seasonal=list(order=c(0, 0, 1), period=7),
include.drift=T,
xreg = as.matrix(data.3[,-c(1,2,3,5)]),
method = "ML"))
#AIC 6853.18
(model3.4 = Arima(data.3$J3,
order=c(1, 1, 1),
seasonal=list(order=c(1, 0, 1), period=7),
include.drift=T,
#xreg = as.matrix(data.3[,-c(1,2,3,5)]),
method = "ML"))
#AIC 6842.05## Series: data.3$J3
## Regression with ARIMA(1,1,1)(0,0,1)[7] errors
##
## Coefficients:
## ar1 ma1 sma1 drift J1L1 J1L2 J2L1 J2L2
## 0.7154 -1.0000 0.0656 0.4232 -0.0304 0.0233 0.1806 -0.1047
## s.e. 0.0298 0.0063 0.0419 0.1056 0.0386 0.0387 0.1183 0.1179
##
## sigma^2 = 8632: log likelihood = -3419.79
## AIC=6857.59 AICc=6857.91 BIC=6896.78
## Series: data.3$J3
## Regression with ARIMA(1,1,1)(0,0,1)[7] errors
##
## Coefficients:
## ar1 ma1 sma1 drift J1L1 J2L1
## 0.7156 -1.0000 0.0671 0.4167 -0.0330 0.1768
## s.e. 0.0298 0.0069 0.0417 0.0969 0.0381 0.1157
##
## sigma^2 = 8615: log likelihood = -3420.22
## AIC=6854.43 AICc=6854.63 BIC=6884.91
## Series: data.3$J3
## Regression with ARIMA(1,1,1)(0,0,1)[7] errors
##
## Coefficients:
## ar1 ma1 sma1 drift xreg
## 0.7133 -1.0000 0.0683 0.3975 0.0932
## s.e. 0.0297 0.0064 0.0417 0.0939 0.0639
##
## sigma^2 = 8611: log likelihood = -3420.59
## AIC=6853.18 AICc=6853.33 BIC=6879.31
## Series: data.3$J3
## ARIMA(1,1,1)(1,0,1)[7] with drift
##
## Coefficients:
## ar1 ma1 sar1 sma1 drift
## 0.7321 -0.9999 0.9917 -0.9711 0.4593
## s.e. 0.0286 0.0037 0.0092 0.0180 0.0991
##
## sigma^2 = 8429: log likelihood = -3415.03
## AIC=6842.05 AICc=6842.2 BIC=6868.18J3 forecasting
model3.forecast1 = forecast(cars.day.fit3, h = 30)
plot(model3.forecast1, xlim = c(75,90))
lines(cars.day$Junction3)
(rmse3.1 = rmse(model3.forecast1$mean, data.3.full[577:606,1]))
# 124.8629
model3.forecast2 = forecast(model3.1, h = 30,
xreg = as.matrix(data.3.full[577:606,-c(1)]))
plot(model3.forecast2, xlim = c(560,610))
lines(data.3.full$J3)
(rmse3.2 = rmse(model3.forecast2$mean, data.3.full[577:606,1]))
# 110.0729
model3.forecast3 = forecast(model3.4, h = 30)
plot(model3.forecast3, xlim = c(560,610))
lines(data.3.full$J3)
(rmse3.2 = rmse(model3.forecast3$mean, data.3.full[577:606,1]))
# 114.9068## [1] 124.8629
## [1] 110.0729
## [1] 114.9068
look at residuals for better fit
#HiTBATS
fit <- tbats(data.1$J1, seasonal.periods = c(7))
fore = forecast(fit, h = 30)
plot(fore, xlim = c(75,90));
lines(ts(data.1.full$J1,frequency = 7))
# cat("model AIC is", fit$AIC)
# AIC 8661.047
rmse.tbats = rmse(data.1.full$J1[577:606], fore$mean)
#148.5043
Log fits
J1 log
(model1.log.1 = auto.arima(log(data.1$J1),
max.p = 5,
max.q = 5,
max.P = 5,
max.Q = 5,
max.order = 10,
max.d = 3,
max.D = 3,
start.p = 1,
start.q = 1,
start.P = 1,
start.Q = 1,
nmodels = 100,
approximation = T))
model1.log.forecast1 = forecast(model1.log.1, h = 30)
plot(model3.forecast3)
#, xlim = c(560,610)
lines(log(data.1.full$J1))
(rmse1.log.1 = rmse(model1.log.forecast1$mean, data.1.full[577:606,1]))
#1790.762## Series: log(data.1$J1)
## ARIMA(5,1,1) with drift
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ma1 drift
## -0.0768 -0.5805 -0.3777 -0.3406 -0.6302 -0.4245 0.0024
## s.e. 0.0371 0.0292 0.0362 0.0295 0.0341 0.0349 0.0010
##
## sigma^2 = 0.01619: log likelihood = 370.64
## AIC=-725.27 AICc=-725.02 BIC=-690.44
## [1] 1790.762
cat("1\n\n\n\n")
auto.arima(cars.day.train$Junction1,
max.p = 5,
max.q = 5,
max.P = 5,
max.Q = 5,
max.order = 10,
max.d = 3,
max.D = 3,
start.p = 1,
start.q = 1,
start.P = 1,
start.Q = 1,
nmodels = 100,
approximation = T)
cat("\n\n2\n\n\n\n")
auto.arima(cars.day.train$Junction1[-c(1,2)],
max.p = 5,
max.q = 5,
max.P = 5,
max.Q = 5,
max.order = 10,
max.d = 3,
max.D = 3,
start.p = 1,
start.q = 1,
start.P = 1,
start.Q = 1,
nmodels = 100,
approximation = T)## 1
##
##
##
## Series: cars.day.train$Junction1
## ARIMA(1,0,1)(0,1,1)[7] with drift
##
## Coefficients:
## ar1 ma1 sma1 drift
## 0.8128 -0.3473 -0.8211 2.1295
## s.e. 0.0386 0.0629 0.0243 0.3114
##
## sigma^2 = 6325: log likelihood = -3310.96
## AIC=6631.91 AICc=6632.02 BIC=6653.65
##
##
## 2
##
##
##
## Series: cars.day.train$Junction1[-c(1, 2)]
## ARIMA(3,1,3)
##
## Coefficients:
## ar1 ar2 ar3 ma1 ma2 ma3
## 0.2874 0.0921 -0.5859 -0.6139 -0.6891 0.6527
## s.e. 0.0636 0.0648 0.0657 0.0503 0.0661 0.0519
##
## sigma^2 = 24274: log likelihood = -3717.2
## AIC=7448.39 AICc=7448.59 BIC=7478.87my.mat1 = matrix(0, nrow = 2, ncol = 3)
colnames(my.mat1) = c("Model 1", "Model 2", "Model 3")
rownames(my.mat1) = c("AIC","RMSE")
my.mat1[1,1] = 6631.91
my.mat1[1,2] = 6604.65
my.mat1[1,3] = 6608.14
my.mat1[2,1] = 114.34
my.mat1[2,2] = 109.64
my.mat1[2,3] = 115.10
knitr::kable(my.mat1)
my.mat2 = matrix(0, nrow = 2, ncol = 3)
colnames(my.mat2) = c("Model 1", "Model 2", "Model 3")
rownames(my.mat2) = c("AIC","RMSE")
my.mat2[1,1] = 5435.15
my.mat2[1,2] = 5415.31
my.mat2[1,3] = 5411.81
my.mat2[2,1] = 43.50
my.mat2[2,2] = 43.36
my.mat2[2,3] = 42.88
knitr::kable(my.mat2)
my.mat3 = matrix(0, nrow = 2, ncol = 3)
colnames(my.mat3) = c("Model 1", "Model 2", "Model 3")
rownames(my.mat3) = c("AIC","RMSE")
my.mat3[1,1] = 6878.07
my.mat3[1,2] = 6857.59
my.mat3[1,3] = 6842.05
my.mat3[2,1] = 124.86
my.mat3[2,2] = 110.07
my.mat3[2,3] = 114.90
knitr::kable(my.mat3)| Model 1 | Model 2 | Model 3 | |
|---|---|---|---|
| AIC | 6631.91 | 6604.65 | 6608.14 |
| RMSE | 114.34 | 109.64 | 115.10 |
| Model 1 | Model 2 | Model 3 | |
|---|---|---|---|
| AIC | 5435.15 | 5415.31 | 5411.81 |
| RMSE | 43.50 | 43.36 | 42.88 |
| Model 1 | Model 2 | Model 3 | |
|---|---|---|---|
| AIC | 6878.07 | 6857.59 | 6842.05 |
| RMSE | 124.86 | 110.07 | 114.90 |
note plot residuals for all models
LateX
Model 1
1
\[ (1-0.81B)(1+0.82B^7)(1-B^7)J_{1,t} = 2.12 + (1-0.34B)w_{1,t} \]
par(mfrow = c(2,2))
plot(cars.day.fit1$residuals, ylab = "residual", main = "residuals")
hist(cars.day.fit1$residuals, freq = F, xlab = "residuals", main = "residual histogram", breaks = 30)
lines(-600:600,dnorm(-600:600, mean = mean(cars.day.fit1$residuals), sd = sd(cars.day.fit1$residuals)), col = "red")
acf(cars.day.fit1$residuals, lag.max = 15, main = "residual ACF")
pacf(cars.day.fit1$residuals, lag.max = 15, main = "residual PACF")
2
\[ (1-0.822B)(1+0.82B^7)(1-B^7)J_{1,t} = 1.906 +(1-0.401B)w_{1,t} + 0.33J_{2,t-1} \]
par(mfrow = c(2,2))
plot(model1.14$residuals, ylab = "residual", main = "residuals")
hist(model1.14$residuals, freq = F, xlab = "residuals", main = "residual histogram", breaks = 30)
lines(-600:600,dnorm(-600:600, mean = mean(model1.14$residuals), sd = sd(model1.14$residuals)), col = "red")
acf(model1.14$residuals, lag.max = 15, main = "residual ACF")
pacf(model1.14$residuals, lag.max = 15, main = "residual PACF")
\[ 1=2-1 \]
Model 2
\[ (1-1.41B+0.43B^2)(1-B^7)J_{2,t} = 0.72 + (1-0.84B)(1-0.79B^7)w_{2,t} \]
par(mfrow = c(2,2))
plot(cars.day.fit2$residuals, ylab = "residual", main = "residuals")
hist(cars.day.fit2$residuals, freq = F, xlab = "residuals", main = "residual histogram", breaks = 30)
lines(-600:600,dnorm(-600:600, mean = mean(cars.day.fit2$residuals), sd = sd(cars.day.fit2$residuals)), col = "red")
acf(cars.day.fit2$residuals, lag.max = 15, main = "residual ACF")
pacf(cars.day.fit2$residuals, lag.max = 15, main = "residual PACF")
\[ (1-1.37B+0.39B^2)(1-B^7)J_{2,t} = 0.64 + (1-0.85B)(1-0.79B^7)w_{2,t} + 0.030J_{1,t-1} + 0.0019J_{1,t-2} + 0.026J_{3,t-1} + 0.0080J_{3,t-2} \]
par(mfrow = c(2,2))
plot(model2.1$residuals, ylab = "residual", main = "residuals")
hist(model2.1$residuals, freq = F, xlab = "residuals", main = "residual histogram", breaks = 30)
lines(-600:600,dnorm(-600:600, mean = mean(model2.1$residuals), sd = sd(model2.1$residuals)), col = "red")
acf(model2.1$residuals, lag.max = 15, main = "residual ACF")
pacf(model2.1$residuals, lag.max = 15, main = "residual PACF")
\[ (1-1.37B+0.39B^2)(1-B^7)J_{2,t} = 0.65 + (1-0.85B)(1-0.79B^7)w_{3,t} +0.029J_{1,t-1} 0.027J_{3,t-1} \]
par(mfrow = c(2,2))
plot(model2.3$residuals, ylab = "residual", main = "residuals")
hist(model2.3$residuals, freq = F, xlab = "residuals", main = "residual histogram", breaks = 30)
lines(-600:600,dnorm(-600:600, mean = mean(model2.3$residuals), sd = sd(model2.3$residuals)), col = "red")
acf(model2.3$residuals, lag.max = 15, main = "residual ACF")
pacf(model2.3$residuals, lag.max = 15, main = "residual PACF")
Model 3
\[ (1-0.69B)\nabla J_{3,t} = (1-0.97B)(1+0.067B^7)w_{3,t} \]
par(mfrow = c(2,2))
plot(cars.day.fit3$residuals, ylab = "residual", main = "residuals")
hist(cars.day.fit3$residuals, freq = F, xlab = "residuals", main = "residual histogram", breaks = 30)
lines(-600:600,dnorm(-600:600, mean = mean(cars.day.fit3$residuals), sd = sd(cars.day.fit3$residuals)), col = "red")
acf(cars.day.fit3$residuals, lag.max = 15, main = "residual ACF")
pacf(cars.day.fit3$residuals, lag.max = 15, main = "residual PACF")
\[ (1-0.71B)\nabla J_{3,t} = 0.42+(1-1B)(1+0.065B^7)w_{3,t} -0.030J_{1,t-1} + 0.023J_{1,t-2} + 0.18J_{2,t-1} -0.10J_{2,t-2} \]
par(mfrow = c(2,2))
plot(model3.1$residuals, ylab = "residual", main = "residuals")
hist(model3.1$residuals, freq = F, xlab = "residuals", main = "residual histogram", breaks = 30)
lines(-600:600,dnorm(-600:600, mean = mean(model3.1$residuals), sd = sd(model3.1$residuals)), col = "red")
acf(model3.1$residuals, lag.max = 15, main = "residual ACF")
pacf(model3.1$residuals, lag.max = 15, main = "residual PACF")
\[ (1-0.73B)(1-0.99B^7)\nabla J_{3,t} = 0.45 + (1-B)(1-0.97B^7)w_{3,t} \]
par(mfrow = c(2,2))
plot(model3.4$residuals, ylab = "residual", main = "residuals")
hist(model3.4$residuals, freq = F, xlab = "residuals", main = "residual histogram", breaks = 30)
lines(-600:600,dnorm(-600:600, mean = mean(model3.4$residuals), sd = sd(model3.4$residuals)), col = "red")
acf(model3.4$residuals, lag.max = 15, main = "residual ACF")
pacf(model3.4$residuals, lag.max = 15, main = "residual PACF")
Relationships between vars in the data.numbers
J1
par(mfrow = c(2,2))
plot(data.1.full$J2L1,data.1.full$J1, main = "J1 | J2(t-1)", xlab = "J2(t-1)", ylab = "J1")
abline(v = 450, col = "red", lty = 2)
plot(data.1.full$J2L2,data.1.full$J1, main = "J1 | J2(t-2)", xlab = "J2(t-2)", ylab = "J1")
abline(v = 450, col = "red", lty = 2)
plot(data.1.full$J3L1,data.1.full$J1, main = "J1 | J3(t-1)", xlab = "J3(t-1)", ylab = "J1")
plot(data.1.full$J3L2,data.1.full$J1, main = "J1 | J3(t-2)", xlab = "J3(t-2)", ylab = "J1")
data.1.full2 = data.1.full
data.1.full2$J2L1L450 = ifelse(data.1.full$J2L1<450,data.1.full$J2L1,0)
data.1.full2$J2L1G450 = ifelse(data.1.full$J2L1>=450,1,0)
#data.1.full2$J2L1L4502 = ifelse(data.1.full$J2L1<450,data.1.full$J2L1,0)
data.1.2 = data.1.full2[-c(577:606),]
(model1.15 = Arima(data.1$J1,
order=c(1, 0, 1),
seasonal=list(order=c(0, 1, 1), period=7),
include.drift=T,
xreg = as.matrix(data.1.2[,-c(1,2,3,4,5)])))
#AIC 6608.14
par(mfrow = c(2,2))
plot(model1.15$residuals, ylab = "residual", main = "residuals")
hist(model1.15$residuals, freq = F, xlab = "residuals", main = "residual histogram", breaks = 30)
lines(-600:600,dnorm(-600:600, mean = mean(model1.15$residuals), sd = sd(model1.15$residuals)), col = "red")
acf(model1.15$residuals, lag.max = 15, main = "residual ACF")
pacf(model1.15$residuals, lag.max = 15, main = "residual PACF")
par(mfrow = c(1,1))
model1.forecast3 = forecast(model1.15, h = 30, xreg = as.matrix(data.1.full2[577:606,-c(1:5)]))
plot(model1.forecast2, xlim = c(560,610))
lines(data.1.full$J1)
rmse1.3 = rmse(model1.forecast3$mean, data.1.full$J1[577:606])
#115.1091## Series: data.1$J1
## Regression with ARIMA(1,0,1)(0,1,1)[7] errors
##
## Coefficients:
## ar1 ma1 sma1 drift J2L1L450 J2L1G450
## 0.8067 -0.3697 -0.8134 1.9971 0.3223 143.3700
## s.e. 0.0413 0.0662 0.0260 0.3110 0.1350 67.5697
##
## sigma^2 = 6300: log likelihood = -3297.07
## AIC=6608.14 AICc=6608.34 BIC=6638.55
\[ (1-0.80B)(1-B^7)J_{1,t} = 1.99+(1-0.36B)(1-0.81B^7)+ 0.32J_{2,t-1}(J_{2,t-1}<450) + 143.37(J_{2,t-1}\ge450) \]
test \(\mathbb{I}\)
\[
\mathbb{I}
\]
J2
par(mfrow = c(2,2))
plot(data.2.full$J1L1,data.2.full$J2, main = "J2 | J1(t-1)", xlab = "J1(t-1)", ylab = "J2")
plot(data.2.full$J1L2,data.2.full$J2, main = "J2 | J1(t-2)", xlab = "J1(t-2)", ylab = "J2")
plot(data.2.full$J3L1,data.2.full$J2, main = "J2 | J3(t-1)", xlab = "J3(t-1)", ylab = "J2")
plot(data.2.full$J3L2,data.2.full$J2, main = "J2 | J3(t-2)", xlab = "J3(t-2)", ylab = "J2")
J3
par(mfrow = c(2,2))
plot(data.3.full$J1L1,data.3.full$J3, main = "J3 | J1(t-1)", xlab = "J1(t-1)", ylab = "J3")
plot(data.3.full$J1L2,data.3.full$J3, main = "J3 | J1(t-2)", xlab = "J1(t-2)", ylab = "J3")
plot(data.3.full$J2L1,data.3.full$J3, main = "J3 | J2(t-1)", xlab = "J2(t-1)", ylab = "J3")
plot(data.3.full$J2L2,data.3.full$J3, main = "J3 | J2(t-2)", xlab = "J2(t-2)", ylab = "J3")
residulal correlations
# J1 model1.14
# J2 model2.3
# J3 model3.4
ggCcf(model1.14$residuals, model2.3$residuals)+ggtitle("CCF Plot for Junctions 1 & 2")
ggCcf(model1.14$residuals, model3.4$residuals)+ggtitle("CCF Plot for Junctions 1 & 3")
ggCcf(model2.3$residuals, model3.4$residuals)+ggtitle("CCF Plot for Junctions 2 & 3")
\[ (1-0.81B)(1-B^7)J_{1,t} = (1-0.34B)(1-0.82B^7)w_{1,t} \]
\[ (1-0.82B)(1-B^7)J_{1,t} = 1.90+ (1-0.40B)(1-0.82B^7)w_{1,t} +0.33J_{2,t-1} \]